/* $Header: mutil.c,v 2.2 01/10/02 10:32:00 przybyls Exp $ */

/***********************************************************************/
/* MODULE NAME:    math                                                */
/* ============                                                        */
/*                                                                     */
/* Author:         Rick                                                */
/*                                                                     */
/* MODULE PURPOSE:  general purpose math routines used for flow/press  */
/*                  calculations (foreground and background)           */
/*                                                                     */
/*                                                                     */
/*                                                                     */
/***********************************************************************/


/*****************************/
/*   GLOBAL INCLUDE FILES    */
/*****************************/

#ifndef dos
#include <m68332.h>
#endif

/*****************************/
/*   LOCAL INCLUDE FILES     */
/*****************************/

#include "typedefs.h"

/*****************************/
/*   GLOBAL VARIABLES        */
/*****************************/


/*****************************/
/*   EXTERNAL VARIABLES      */
/*****************************/


/*****************************/
/*   STATIC VARIABLES        */
/*****************************/
/**************************************************/
/* the following table was created with mk_sqrt.c */
/**************************************************/

const U_INT16 sqrt_tbl[] = {

/*0x1000000 */   4096,  5792,  7094,  8192,  9158, 10034, 10836, 11586,
/*0x9000000 */  12288, 12952, 13584, 14188, 14768, 15326, 15864, 16384,
/*0x11000000*/  16888, 17378, 17854, 18318, 18770, 19212, 19644, 20066,
/*0x19000000*/  20480, 20886, 21284, 21674, 22058, 22434, 22806, 23170,
/*0x21000000*/  23530, 23884, 24232, 24576, 24914, 25250, 25580, 25906,
/*0x29000000*/  26228, 26546, 26860, 27170, 27476, 27780, 28080, 28378,
/*0x31000000*/  28672, 28964, 29252, 29536, 29820, 30100, 30376, 30652,
/*0x39000000*/  30924, 31194, 31462, 31728, 31990, 32252, 32510, 32768,
/*0x41000000*/  33024, 33276, 33528, 33776, 34024, 34270, 34514, 34756,
/*0x49000000*/  34996, 35236, 35472, 35708, 35942, 36174, 36406, 36636,
/*0x51000000*/  36864, 37090, 37316, 37540, 37764, 37984, 38204, 38424,
/*0x59000000*/  38642, 38858, 39074, 39288, 39500, 39712, 39922, 40132,
/*0x61000000*/  40340, 40548, 40754, 40960, 41164, 41368, 41570, 41772,
/*0x69000000*/  41972, 42170, 42370, 42566, 42764, 42960, 43154, 43348,
/*0x71000000*/  43542, 43734, 43924, 44116, 44306, 44494, 44682, 44870,
/*0x79000000*/  45056, 45242, 45426, 45612, 45794, 45978, 46160, 46340,
/*0x81000000*/  46522, 46702, 46880, 47060, 47238, 47414, 47592, 47768,
/*0x89000000*/  47942, 48118, 48292, 48464, 48638, 48810, 48982, 49152,
/*0x91000000*/  49322, 49492, 49662, 49830, 49998, 50166, 50332, 50498,
/*0x99000000*/  50664, 50830, 50994, 51160, 51322, 51486, 51648, 51810,
/*0xa1000000*/  51972, 52134, 52294, 52454, 52614, 52774, 52932, 53090,
/*0xa9000000*/  53248, 53406, 53562, 53718, 53874, 54030, 54184, 54340,
/*0xb1000000*/  54494, 54648, 54800, 54954, 55106, 55258, 55410, 55560,
/*0xb9000000*/  55712, 55862, 56012, 56162, 56310, 56460, 56608, 56756,
/*0xc1000000*/  56904, 57050, 57198, 57344, 57490, 57636, 57782, 57926,
/*0xc9000000*/  58070, 58216, 58360, 58502, 58646, 58788, 58932, 59074,
/*0xd1000000*/  59216, 59356, 59498, 59638, 59780, 59920, 60060, 60198,
/*0xd9000000*/  60338, 60476, 60616, 60754, 60892, 61030, 61166, 61304,
/*0xe1000000*/  61440, 61576, 61712, 61848, 61984, 62118, 62254, 62388,
/*0xe9000000*/  62522, 62656, 62790, 62924, 63058, 63190, 63322, 63454,
/*0xf1000000*/  63588, 63718, 63850, 63982, 64112, 64244, 64374, 64504,
/*0xf9000000*/  64634, 64764, 64892, 65022, 65150, 65280, 65408, 65535 };



   const  U_INT16 point_seven_power_table[] = {
         0,     48,     78,    104,    127,    149,    170,    189,
       207,    225,    243,    259,    276,    292,    307,    322,
       337,    352,    366,    380,    394,    408,    422,    435,
       448,    461,    474,    487,    499,    512,    524,    536,
       548,    560,    572,    584,    595,    607,    618,    630,
       641,    652,    663,    674,    685,    696,    707,    718,
       728,    739,    749,    760,    770,    781,    791,    801,
       811,    822,    832,    842,    852,    861,    871,    881,
       891,    901,    910,    920,    930,    939,    949,    958,
       968,    977,    986,    996,   1005,   1014,   1023,   1033,
      1042,   1051,   1060,   1069,   1078,   1087,   1096,   1105,
      1114,   1122,   1131,   1140,   1149,   1157,   1166,   1175,
      1184,   1192,   1201,   1209,   1218,   1226,   1235,   1243,
      1252,   1260,   1269,   1277,   1285,   1294,   1302,   1310,
      1318,   1327,   1335,   1343,   1351,   1359,   1367,   1376,
      1384,   1392,   1400,   1408,   1416,   1424,   1432,   1440,
      1448,   1456,   1463,   1471,   1479,   1487,   1495,   1503,
      1510,   1518,   1526,   1534,   1541,   1549,   1557,   1564,
      1572,   1580,   1587,   1595,   1603,   1610,   1618,   1625,
      1633,   1640,   1648,   1655,   1663,   1670,   1678,   1685,
      1692,   1700,   1707,   1715,   1722,   1729,   1737,   1744,
      1751,   1759,   1766,   1773,   1780,   1788,   1795,   1802,
      1809,   1816,   1824,   1831,   1838,   1845,   1852,   1859,
      1866,   1874,   1881,   1888,   1895,   1902,   1909,   1916,
      1923,   1930,   1937,   1944,   1951,   1958,   1965,   1972,
      1979,   1986,   1993,   1999,   2006,   2013,   2020,   2027,
      2034,   2041,   2047,   2054,   2061,   2068,   2075,   2081,
      2088,   2095,   2102,   2109,   2115,   2122,   2129,   2135,
      2142,   2149,   2155,   2162,   2169,   2175,   2182,   2189,
      2195,   2202,   2209,   2215,   2222,   2228,   2235,   2242,
      2248,   2255,   2261,   2268,   2274,   2281,   2287,   2294,
      2300,   2307,   2313,   2320,   2326,   2333,   2339,   2346 };

/*****************************/
/*   IMPLEMENT SECTION       */
/*****************************/

/*****************************************************************************/
/* FUNCTION:       Sqrt_LU( n )                                              */
/* PURPOSE:        lookup the square root of a 32bit unsigned integer        */
/* RETURNED VALUE: a 16bit unsigned integer approximating the sq root of n.  */
/* LIMITATIONS:    Errors in this routine have been numerically calculated   */
/*                 to be less than +5/-2 relative to the newtons method for  */
/*                 32bit numbers and less than +2/-1 for 16bit numbers.      */
/* GLOBALS USED:   sqrt_tbl  defined in this file                            */
/* CALLS:          tableU  defined locally and by <m68332.h>                 */
/*                                                                           */
/* First, the format of the lookup table, sqrt_tbl.  Entries represent the   */
/* square roots of the index+1 multipled by 0x01000000 (4096 squared).       */
/*                                                                           */
/*     entry 0     ( 0x01000000 ) ^.5    (0x0100)^.5 * 256    4096           */
/*     entry 1     ( 0x02000000 ) ^.5    (0x0200)^.5 * 256    5792           */
/*     entry 2     ( 0x03000000 ) ^.5    (0x0300)^.5 * 256    7094           */
/*      ...             ...                                     .            */
/*     entry y     ( (y+1)*4096^2 ) ^.5 (y+1*100)^.5 * 256      .            */
/*      ...             ...                                     .            */
/*     entry 253   ( 0xfe000000 ) ^.5    (0xfe00)^.5 * 256   65280           */
/*     entry 254   ( 0xff000000 ) ^.5    (0xff00)^.5 * 256   65408           */
/*     entry 255   (0x100000000 ) ^.5   (0x10000)^.5 * 256   65535           */
/* the last entry (255) is a 16bit aproximation for 0x100000000^.5           */
/*                                                                           */
/* To find the sqrt of n, shift n (multiply or divide) by a constant Z to    */
/* place the most significant 8 bits of the number n in the index field      */
/* with the next most significant 8 bits in the round field.  Since the      */
/* table uses entry number y as the sq root of y+1, subtract one from the    */
/* index.  Use the tblU instruction to perform the lookup and interpolate,   */
/* then shift (divide or multiply) by the sq root of the constant Z.  Then   */
/* divide by 256 to correct for the table scaling.                           */
/*                                                                           */
/* The reason for shifting the entries in the table (y -> y+1) is that       */
/* the algorithm used by the tblU instruction looks at entry and entry+1     */
/* to perform the interpolation.  Shifting the entries avoids going          */
/* beyond the last entry when doing interpolations between 0xff000000        */
/* and 0xffffffff.  This entry shifting thus requires a test for the         */
/* square root of zero.                                                      */
/*                                                                           */
/* The reason for shifting by Z is to get the most significant bits in       */
/* the index field to get better resolution for smaller numbers.             */
/* Algebraically, this works out to be just multiplying and dividing by      */
/* constants.                                                                */
/*                                                                           */
/*          shift by a       decrement  do the  undo the    undo table       */
/*          constant Z       the index  lookup   shift Z    scaling (256)    */
/*  00000009   << 12   9000    8f00     c000     >> 6   0300   >> 8    0003  */
/*  000000ff   <<  8   ff00    fe00     ff7f     >> 4   0ff7   >> 8    000f  */
/*  0000083f   <<  4   83f0    82f0     b7c7     >> 2   2df1   >> 8    002d  */
/*      n     / * Z    n*Z      (n*Z)^.5*256  / * Z^.5  n^.5*256       n^.5  */
/*  0000cbfe           cbfe    cafe     e485            e485   >> 8    00e4  */
/*  000f1fff   >>  4   f1ff    f0ff     f8e6     << 2  3e389   >> 8    03e3  */
/*  00fade21   >>  8   fade    f9de     fd69     << 4  fd690   >> 8    0fd6  */
/*  0c0de000   >> 12   c0de    bfde     de33     << 6 378cc0   >> 8    378c  */
/*  deadbeaf   >> 16   dead    ddad     eec2     << 8 eec200   >> 8    eec2  */
/*****************************************************************************/

U_INT16 SqrtLU( U_INT32 num_sqrd )

{

#define DEC_INDEX(num) (((((num)&0xff00)-1)&0xff00)|((num)&0x00ff))


     if ( num_sqrd <= 0x0000ffff )
     {
          if ( num_sqrd <= 0x000000ff )
          {
               if ( num_sqrd <= 0x0000000f )
               {
                   if ( num_sqrd == 0 )
                   {
                        return ( 0 );
                   }
                   else
                   {
                        return(tableU((U_INT16 *) sqrt_tbl, DEC_INDEX(num_sqrd<<12))>>14);
                   }
               }
               else
               {
                        return(tableU((U_INT16 *)sqrt_tbl, DEC_INDEX(num_sqrd<< 8))>>12);
               }
          }
          else
          {
               if ( num_sqrd <= 0x00000fff )
               {
                        return(tableU((U_INT16 *)sqrt_tbl, DEC_INDEX(num_sqrd<< 4))>>10);
               }
               else
               {
                        return(tableU((U_INT16 *)sqrt_tbl, DEC_INDEX(num_sqrd    ))>> 8);
               }
          }
     }
     else
     {
          if ( num_sqrd <= 0x00ffffff )
          {
               if ( num_sqrd <= 0x000fffff )
               {
                        return(tableU((U_INT16 *)sqrt_tbl, DEC_INDEX(num_sqrd>> 4))>> 6);
               }
               else
               {
                        return(tableU((U_INT16 *)sqrt_tbl, DEC_INDEX(num_sqrd>> 8))>> 4);
               }
          }
          else
          {
               if ( num_sqrd <= 0x0fffffff )
               {
                        return(tableU((U_INT16 *)sqrt_tbl, DEC_INDEX(num_sqrd>>12))>> 2);
               }
               else
               {
                        return(tableU((U_INT16 *)sqrt_tbl, DEC_INDEX(num_sqrd>>16)));
               }
          }
     }

}

/*****************************************************************************/
/* FUNCTION:       Point_Seven_Power(n)                                      */
/* PURPOSE:        calculate n raised to the .7 power                        */
/* RETURNED VALUE: a 16bit unsigned integer approximating n to the .7        */
/* LIMITATIONS:    n <= 0xff00                                               */
/* GLOBALS USED:                                                             */
/* CALLS:                                                                    */
/*****************************************************************************/

U_INT16 PointSevenPower( U_INT16 number )
{
      return tableU((U_INT16 *)point_seven_power_table, number );

}


